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ABSTRACT 

We present results from three-dimensional, self-gravitating, radiation-hydrodynamic simulations of 
low-mass protostellar outflows. We construct synthetic observations in ^^CO in order to compare 
with observed outflows and evaluate the effects of beam resolution and outflow orientation on inferred 
outflow properties. To facilitate the comparison, we develop a quantitative prescription for measur- 
ing outflow opening angles. Using this prescription, we demonstrate that, in both simulations and 
synthetic observations, outflow opening angles broaden with time similarly to observed outflows. How- 
ever, the interaction between the outflowing gas and the turbulent core envelope produces significant 
asymmetry between the red and blue shifted outflow lobes. We find that applying a velocity cutoff 
may result in outflow masses that are underestimated by a factor 5 or more, and masses derived from 
optically thick CO emission further underpredict the mass of the high-velocity gas by a factor of 5-10. 
Derived excitation temperatures indicate that outflowing gas is hotter than the ambient gas with 
temperature rising over time, which is in agreement with the simulation gas temperatures. However, 
excitation temperatures are otherwise not well correlated with the actual gas temperature. 
Subject headings: stars: formation stars: outflows 



1. INTRODUCTION 

Young protostars power high- velocity jets that en- 
train and unbind a large fraction of the natal proto- 
stellar core gas. Outflow mass rates are estimated to 
be co niparable in ma gs:nitude to protostellar accretion 
rates (jBontemps et al. 1996). Consequently, outflows 
likely play an important role in removing core mass and 
. term i nating; the a ccretion process ( Matzner fc McKe^ 
■ I2QQQI : lMversll2QQ9h . The largest outflows accelerate gas 
to velocities exceeding 100 kms~^ and extend across 
several parsecs. On these scales, outflows powerfully 
impact their environment and potentially inject sig- 
nificant energy bac k into the parent molecular cloud 
(|Matzner fc McKeel[r999l This feedback may be par- 
tially respo nsible for maintaining; turbulence on 0.1- 
1 PC scales (Nakamura & Li" '2007^ 'Sw ift fc Weld] [2008[ 
lArce fc Sargent 2005; Arce et al. 2010). 

Due to observational resolution limits and the inher- 
ently opaque nature of young protostellar cores, the 
central outflow engine is not well understood. Out- 
flow morphology, velocity distribution, and opening an- 
gle provide important constraints for theoretical mod- 
els. For example, observations suggest that outflows 
develop fro m a narrow jet- like niQrpholog;y into a wide - 
angle wind (|Arce fc Sargentll2006l : ISeale fc Loo"iievll2008f ) . 
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Outflow velocities also follow a Hubble-like relationship, 
with the gas vel ocity increasing linearly from the source 
(|Arce et al.l 120071) . Among proposed models, the wide- 
angle- windmodey^i fc Shu 1996) and the jet bow-shock 
model (|Chernin fc Massonl ll995) ap pear to best exp lain 
the observed outflow characteristics (|Lee et al.ll2000f ). 

However, a number of uncertainties complicate inter- 
pretation of the observations. Orientation of the disk- 
outflow geometry relative to the line-of-sight mak e pro- 
tostell ar age es timations imprecise (|Ladd et al.l 119981 
.Robitaille et al.|[2QQ6.) . Once the outflow gas shocks and 
cools, it quickly mixes with ambient turbulent gas and 
becomes impossible to distinguish. This means that ob- 
served outflow gas may only span the previous few thou- 
sand years of the flow. In addition, the ambient low- 
density, high-linewidth cloud gas {c^ 10 kms~^) along 
the line-of-sight makes identification of a wide -angle, low- 
veloc ity outflow component problematic (e.g.. I Arce et al.l 

Eoia). 

Precession of the jet due to rotational wobbling of the 
protostar m ay increase the broadn ess and dumpiness of 
the outflow (|R.osen fc Smithll2004f ). It is likely that vari- 
ability in the protostellar a ccretion rate causes variability 
in outflow properties (e .g., ICernicharo fc Reipurt hi [19961 : 
I Arce fc Goodmanll2002l ). Consequently, bursty accretion 
may produce a signature in the velocity spectrum and 
outflow morphology. However, this is difficult to demon- 
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strate observationally. 

In this work, we use gravito-radiation-hydrodynamic 
simulations modeling protostellar outflows to study out- 
flow evolution as a function of time. First, we intro- 
duce a quantitative method for determing opening an- 
gles and use this to characterize the outflow proper- 
ties of the raw simulation data. Next, we investigate 
the dependence on inclination and observed resolution. 
Finally, we present synthetic observations of the simu- 
lations in ^^CO. Throughout, we compare with inter- 
ferometric observations of s even protostellar outflows 
studied by lArce fc Sargentl (j2QQ6f ) (henceforth AS06). 
These observations offer a high-resolution picture of the 
inner outflow regions, the interaction between outflow 
and envelope, and the evolution of outflow characteris- 
tics over time. Since the simulations provide complete 
three-dimensional information, they allow us to compare 
position-position-position (ppp) data with synthetic and 
observational position-position-velocity (ppv) data and 
examine the observations. 

We outline the simulation methodology and initial con- 
ditions in §2. We define our procedure for characterizing 
the outflows in §3. We present the data analysis, syn- 
thetic observations and observational comparison in §4. 
We summarize our results in §5. 

2. NUMERICAL METHODS 

The simulations are performed using the ORION 
Adaptive Mesh Refinement (AMR) code. For our 
study, we use the ra diation- hydro d ynamic s, self- 
gravitating simulations of lOffner et al.l ()2QQ9bl ). hence- 
forth OKMK09, as initial conditions. 

The OKMK09 simulations have a mean density of 
4.46 X 10~^^ g cm-^, 3D Mach number of 6.6, and a total 
mass of 185 M©. Energy is injected at a constant rate to 
maintain the level of turbulence (e.g. Stone et al. 1998). 
Particles are inserte d in regions of the flow that exceed 
the Jeans condition (|Krumholz et al.ir2QQ4f ). These stars 
are endowed with a sub-grid stellar evolution model that 
includes accretion luminosity down to the stellar surface, 
Kelvin-Helmholz contraction, and nuclear burning. The 
original calculations include the effects of radiation feed- 
back from the forming stars but neglect protostellar out- 
flows, which is the subject of our investigation. Note that 
we use the terms outflow and wind interchangeably and 
do not distinguish between gas that is directly launched 
near the protostar and gas that is swept up and entrained 
by this gas. 

From the OKMK09 simulations, we select two form- 
ing protostars to receive additional refinement beginning 
prior to their formation. In these high-resolution regions, 
the minimum cell size is 4AU. The remainder of the sim- 
ulation is evolved with the previous number of 4 AMR 
levels, i.e., 32 AU minimum cell size, and without out- 
flows. Although more computationally expensive than 
modeling isolated protostars, this method allows us to 
use self-consistent turbulent initial conditions and follow 
the evolution of the outflows within a turbulent, clus- 
tered star- forming environment. Henceforth, we refer to 
the two high-resolution runs as Rl and R2. 

To the selected protostars, we add a sub-grid model 
for pr otostellar winds bas ed up on Matzner & McKee 
unningham et al.l (l2QllD describe the details of 
the model implementation in ORION, which we briefly 
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Fig. 1. — Distribution of injected angular momentum, ^(6,6o) 
versus polar angle, 6, for = 0.1. 

summarize here. The outflow model is characterized by 
three dimensionless parameters that specify the outflow 
ejection efficiency, outflow velocity, and momentum dis- 
tribution. The mass ejection rate, fw, gives the frac- 
tion of infalling gas that is accelerated into a wind. 
This f raction is observationally uncertain, bu t the disk 
wind (jPelletier fc Pudritzl IT992t ) and X-wind (|Shu et al.l 
Il994f ) models predict 0.1 — 0.33. Here, we adopt 

fw = 0.33. Consequently, 1.0/(1 + f^) of the infalling 
gas accretes onto the star, while /^/(l + /w) is launched 
in an outflow. 

The wind launching velocity is given by the Keplerian 
velocity at the stellar surface, vk = yGM^Ji\. In the 
case of high-mass protostars these velocities can exceed 
200 kms~^, greatly constraining the numeric al timestep 
of the calculation. ICunningham et ahl (j2Qllf ) solve this 
problem by limiting the outflow velocity to a fraction of 
the Keplerian speed, In the calculations we present 
here, we set fv = ^ since the stars forming are low- mass 
with Keplerian velocities < 100km s~^. 

The direction of the outflow ejection is set by the direc- 
tion of the angular momentum vector of the protostar. 
This is determined by the angular momentum of the ac- 
creting gas, a quantity that depends upon the turbulent 
properties of the core and evolves over the calculation. 
Thus, the direction of the wind is not fixed or set as an 
input parameter but is self-consistently dictated by the 
hydrodynamic evolution of the accretion flow. As a con- 
sequence, the launching axis is roughly parallel to the 
angular momentum vector of the accretion disk. 

We def ine the effective ope ning angle of the wind, Oq, 
following iMatzner fc McKed (|l999). They characterize 
the angular distribution of the outflow momentum: 



In 



(1) 



where is the polar angle measured from the proto- 
star's rotation axis. Through comparisons wit h obser- 
vation al data of low-mass protostars Ma tzner fc McKed 
(|l999f ) flnd that ^ 0.05 and suggest a flducial value of 
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^0 = 0.01. This value results in a strongly peaked dis- 
tribution, which numerically deposits nearly all the mo- 
mentum in a couple cells along the rotational axes. Here, 
we adopt a slightly larger value of = 0.1 as shown in 
Figure [U which is comparable to the angular width of 
an jndividual fine cell. The momentum assigned to a cell 
is ^, the momentum distribution function averaged over 
the polar angle subtended by that cell. 

The wind is injected into cells with radial distances 
4Ax < r < SAx from the protostar, so that the wind in- 
jection zone lies outside the accretion region. This serves 
the dual purpose of providing better resolution in the in- 
jection region and allowing gas to continue to accrete 
onto the central protostar. In practice, <f is also set to 
when becomes close to 7r/2. Our outflow algorithm is 
fully mass conserving. Gas that was previously automat- 
ically accreted onto the stars in the original simulation 
is instead divided between accreted gas and outflow gas 
that is deposite d back onto the num erical grid. 

Instead of the Pollack et al.' (1994) dust opacity model 
used by 0KMKQ9, w e switch to an updated model from 
[Semenov et al.l (|2003f ). which assumes a standard iron 
abundance and treats the grains as composite aggregates. 
The strong bow shocks produced by outflowing gas run- 
ning into ambient material can generate temperatures 
well in excess of the dust destruction temperature. Since 
the outflows simulated here are young and very low-mass, 
only a small number of cells ever reach temperatures 
above 1000 K. Nonetheless, we inc lude the treatment 
for atomic line cooling described by "Cunningham et alj 
(2011), which implicitly solves for cell temperatures ex- 
ceeding 10^ K. 

3. ANALYSIS 
3.1. Outflow Identification 

Identifying the gas associated with outflowing mate- 
rial remains a key observational challenge. Generally, 
outflow gas is higher velocity, hotter, and less dense than 
the accreting material. It is difficult to use temperature 
as a means to distinguish outflow gas. In nearby well- 
resolved regions, it is possible to image outflow cavities 
in scattered light (|Seale fc Loonevir2008f ). However, ob- 
servers typically use low-density tracers such as ^^CO to 
map the outflow gas, since this data supplies both di- 
rect velocity information and indirect estimates of the 
gas mass. Within molecular line data cubes, line-of-sight 
velocity remains the primary means of identification. In 
AS06, outflows are defined by selectively integrating over 
the CO (1-0) emission in a range of high- velocity chan- 
nels determined by visual inspection. 

To identify the outflows in the simulations, we set a 
minimum outflow gas velocity of 2 kms~^. We then gen- 
erate a column density map of cells above this cutoff, 
which is analogous to an observed intensity map inte- 
grated over selected velocity channels. 

Figure [2] shows the angle distributions of integrated 
emission for cartoon elliptical outflows. In this 2D space, 
we treat each outflow lobe independently and character- 
ize outflow morphology using five parameters: opening 
angle, length, width, shape, and inclination. To measure 
the opening angle, we calculate the angle with respect 
to the z axis for each pixel that is above a given veloc- 
ity minimum. In the analysis, we weight all such pixels 
equally. The opening angle is then defined to be the full- 
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Fig. 2. — Left: Idealized elliptical outflow shape with half angle 
atan(6/a) = 20 degrees, where a and b are the length and width 
of the major and minor axes, progressively truncated towards the 
center. The black lines show the measured opening angle. Right: 
Angle distribution of the pixels constituting the outflow, where the 
vertical lines indicate the opening angle deflned as the full-width 
quarter maximum of the distribution. 

width-quarter-maximum of the angle distribution. This 
definition is empirically selected to correspond to opening 
angles derived from intensity maps by eye (and protrac- 
tor). We derive the outflow inclination with respect to 
the vertical axis by fitting the angle distribution with ei- 
ther a Gaussian or higher order polynomial. Note that 
an outflow may also be inclined along the line of sight to- 
wards (or away from) the observer. However, our method 
only fits for the inclination in the plane of the sky. For 
an elliptical outflow, the inclination is the location of 
the distribution maximum, while for a more conically 
shaped outflow, it is the local minimum. Note that since 
the simulation inclination is obtained from the projected 
distribution of cells it may not identically correspond to 
the direction in which the outflow is launched. In addi- 
tion, interaction between the outflow and the envelope 
may affect the inferred outflow inclination. 

For a given distribution of connected high-velocity 
cells, the physical length of the outflow is the longest 
extent measured from the center outwards along the out- 
flow major axis. The width is the maximum extent along 
the axis perpendicular to the outflow length. The Gaus- 
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sian parameter, Gauss, is the value of a Gaussian fit 
to the angle distribution. As illustrated by Figure [21 the 
distribution of pixels for elliptically shaped outflows is 
well described by a Gaussian. Our opening angle def- 
inition is formulated to be consistent with an elliptical 
outflow such that atan( Width/Length) 6. As we show 
in ^4.2[ the pixels close to the protostar are likely to be 
dominated by the beam characteristics, and thus are not 
suitable for an opening angle determination. 

3.2. Opening Angle Measurement 

We first characterize an idealized, symmetric outflow 
geometry. Figure [2] shows an ellipse of length 0.1 pc, 
which has been superimposed upon the AMR grid cell 
hierarchy of a simulation output. The effect of truncat- 
ing the outflow length is apparent in the angle distribu- 
tion evolution shown on the right. The inferred angle 
increases with increasing truncation as the wide angle 
cells at the base of the ellipse dominate the distribution. 

Figure [3] shows the opening angle and major axis incli- 
nation for a series of truncated elliptical outflows. The er- 
ror bars indicate the y/N error, where is the number of 
beams contained in the ellipse assuming that the outflow 
is placed 250 pc away and observed with a beam of 4" . 
Figure [3^ illustrates that the error bars do not necessarily 
contain the inferred opening angle for the outflow with 
full information (i.e., no truncation) even though the un- 
certainty increases as the outflow length shrinks. This 
indicates that a poorly resolved outflow may be mea- 
sured to have a fundamentally different opening angle 
than the same outflow with complete information. The 
outflow major axis remains fixed during the truncation. 
Figure [SJd demonstrates that the derived direction of the 
outflow axis is fairly insensitive to the truncation. 

For comparison, we re-analyze the AS06 data to ob- 
tain the opening angles quantitatively. Figure H] shows 
the measure angles for the upper and lower lobes and 
the mean of the two. The slope of the fit we find to 
the data, 0.15, is quite close to the value of 0.16 ± 0.4 
found by AS06. The intercept of 1.22 is also within er- 
ror of the previously published value of 1.1 ±0.2. This is 
encouraging since it suggests that our method is a good 
quantitative alternative to fitting the outflow angles by 
eye. One caveat to the fitting method is that for small 
numbers of pixels the algorithm becomes sensitive to the 
bin size. We indicate this uncertainty using vertical error 
bars on the plot, which are twice the bin size used in the 
fit. 

As shown in the bottom panel of Figure IH the an- 
gle distribution tends to evolve from an elliptical shape 
(Gauss (x^) < 20) to a more conical or irregular shape 
(Gauss (x^) > 20). This trend may be partially an ar- 
tifact of how the outflow is observationally sampled at 
late times rather than an indication of how the geometry 
actually changes. This is supported by the apparent de- 
crease of outflow length with time, which is most likely 
because the older gas has blended with the ambient gas. 

4. RESULTS 

4.1. Outflow Morphology 

Figure [5] shows a volume rending of the simulation ve- 
locities for both runs. The outflow axis in Rl is almost 
directly aligned with the z axis so that Vz traces the 
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Fig. 3. — Top (a): Opening angle as a function of ellipse 
length included in the angle determination. (Ellipses with frac- 
tional lengths 0.2, 0.4, 0.6, 0.8 and 1.0 are shown in Figure O. 
Bottom (b): Inclination of the ellipse major axis relative to the y 
axis as a function of the ellipse length fraction. The gray shading 
shows the range for the non-truncated case. 

fastest moving gas, which ranges up to 70 kms-^ R2 
forms two protostars with fairly distinct outflows. The 
secondary protostar arises from fragmentation in the core 
rather than disk fragmentation, so that the outflow axes 
of the two are somewhat misaligned. However, the two 
R2 protostars are sufficiently close that separating the 
outflow lobes in projection is difficult. Consequently, in 
the following analysis we will use Rl to explore outflow 
evolution for individual sources, while R2 will illustrate a 
case in which observations are confused by the presence 
of a second, unidentifled outflow source. 

All three outflows shown in Figure [5] exhibit a sig- 
niflcant amount of asymmetry between the lobes. The 
sub-grid model launches the outflow gas from the grid 
symmetrically about the protostars. (On the scale of 
the flgure the outflow launching region is contained in 
the central pixel.) This means that any asymmetry that 
arises is directly the result of the interaction between the 
outflow and the turbulent envelope. For example, the 
Rl protostar forms in a more fllamentary structure and 
the net angular momentum direction is almost directly 
aligned with the fllament axes. As a result, the positive 
z outflow lobe is diverted from the launching axis and 
conflned by the dense fllament gas. The outflow lobe 
along the negative z axis breaks out of the fllament and 
extends further. In R2 asymmetry in the core also leads 
to mismatched outflow lobe sizes. 
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Fig. 4. — Opening angle (top) and Gauss parameter (bottom) 
for the 'Arce & Sargent (2006) sample as a function of estimated 
age. Stars indicate the original values derived by Arce &; Sargen^ 
f2006) with the fit shown by the dot-dashed line. The vertical error 
bars on the average opening angle are twice the angle bin size used 
for the fit; in some cases they are smaller than the symbol. The 
dashed line shows the revised fit of the mean opening angle derived 
for our angle definition. 

4.2. Opening Angle Evolution 

Figure [6] shows the Rl outflow evolution as a function 
of time. The outflow gas is identified using a velocity cut- 
off of 2 kms~^, where the full 3D information is used to 
determine the outflow and envelope masses. These cells 
are projected along the x-direction so that the outflow 
opening angle, width, and length is calculated when the 
outflow axis is nearly parallel to the z-axis. A fit to the 
opening angles demonstrates that the lobes widen with 
time similarly to observations. The earliest time only 
contains a small number of cells in the outflow so that 
it is artificially broadened, but a strong trend is appar- 
ent for 3.7 < log(t) < 4.7. (See Section 5 for discussion 
of the origin of this broadening.) However, the differ- 
ent outflow lobes in Figure [6] differ from one another and 
display significant shorter-time variation. For example, 
the upper outflow lobe growth stalls at late times. This 
suggests that even if outflows generally evolve with time, 
variation in individual outflow behavior can be large. 

Considering the large uncertainties in observed out- 
flow ages, poorly constrained line-of-sight inclinations, 
and instrument limitations, the agreement between the 
simulations and the observations shown in Figure |4] is 



Fig. 5. — Volume rendering of the outfiow gas velocity (km s"-*^) 
for Rl at 40 kyr (top) and for R2 at 20 kyr (bottom). The box is 
0.15 pc on a side with the stellar positions marked by green crosses. 

striking. This supports the finding of AS06 that out- 
flow angles evolve with time. However, the protostellar 
masses in these calculations are ^0.1 for the time 
of comparison, so they are still very young. Outfiow 
sizes are < O.lpc, although this is similar to the sizes 
of those in AS06, which are ~ 0.005 — O.lpc. Figure [6] 
illustrates the length of only the connected cells with ve- 
locities > 2km s~^. These cells are also very low-density 
and have a total mass of ~ 0.001 M©. 

The Rl opening angle extrapolated at t = is lower 
than that found by AS06, although the positive z fit 
is still within the observational error. Smearing due to 
the observational beam likely accounts for some of this 
discrepancy. Figure [3 shows the measured opening an- 
gles for Rl as a function of convolution with different 
beam sizes. For a 4" beam, which is comparable to the 
resolution of AS06, some angles may be broadened by 
up to ~80%. Broadening becomes more significant for 
larger beam sizes. The earliest times in the simulation, 
for which the outflow is already broadened by the small 
number of cells, appear least affected by the beam smear- 
ing. 
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Fig. 6. — Opening angles, outflow length, outflow width, and mass for the simulation as a function of time, where the angle are calculated 
for the integrated column of the cells belonging to the outflow as deflned by a velocity > 2km s"-*^. The simulation data has been divided 
into 45 logarithmically spaced bins. The vertical error bars indicate the error given the number of beams contained in the outflow if it is 
250 pc away and observed with a 4" beam. The core mass is deflned as the sum over cells with densities > 10^ g cm~^. 



4.3. Synthetic Observations 

While comparison between the simulation and ob- 
served ^^CO data in ^4.21 is suggestive, it is important 
to compare the data directly with synthetic observations 
in ^^CO. In t his section we use t he radiative transfer code 
MOLLIE (K eto fc Rvbickil[2Ql Q) to model the emission 
from the first seven ^^CO and ^ CO transitions. We as- 
sume a standard abundance of 5.6 x 10~^ ^^CO relative 
to H2 and an abundance of 7.3 x 10~^ ^^CO relative to 
H2, which are comparable to CO es timated abundances 
in low -mass envelopes and outflows (|Carolan et al.|[2QQ8l 
l2QQ9f ). However, there are a wide range of inferred CO 
abundances among different cores and star-forming re- 
gions, values of which are often individually uncertain 
by a factor of two or more. These uncertainties com- 
plicate the estimation of gas mass beyond the already 
challenging problem of identifying the outflow gas and 
accounting for projection. 

We perform the radiative transfer calculation on a cube 
of side 0.15 pc centered on the protostar. We regrid the 
AMR data to produce a grid of 128^ with 130 AU ceh 



resolution nested within a grid of 128^ and 260 AU cell 
resolution. 

Figure [8] shows integrated CO (1-0) emission maps 
from observations of Rl at different inclinations with re- 
spect to the line-of-sight. This includes all gas with ve- 
locities in the range ±20 kms~^, which encompasses the 
minimum and maximum outflow velocity at 130 AU res- 
olution.^ The outflow structure is clearly visible in the 
integrated maps. Such structure may not be apparent 
in observed integrated maps due to emission from ambi- 
ent cloud gas between the source and the observer. The 
right panels show the map an observer would make when 
integrating over gas in higher velocity channels {\v\ > 2 
km s""*^). 

The lower outflow lobe shows substructure that could 
be interpreted as the result of accretion bursts. How- 
ever, the accretion history, and hence the outflow mass 

^ In the central launching region in which the AMR resolution 
is 4 AU, velocities reach ~ 100 kms"-*^. These velocities occur in 
a small volume and are reduced by the flattening to constant pixel 
size. 
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Fig. 7. — Opening angle versus beamsize assuming that the 
outflow is observed at a distance of 250 pc. The color scale indi- 
cates the age of the outflow in units of 10^ years. Symbols for the 
different lobes have been offset for clarity. 

ejection rate, is fairly smooth (see Figure [12]). Instead, 
the discontinuous morphology results from the outflow 
interacting with the turbulent and asymmetric core gas. 

Figure E] shows the opening angle versus inclination 
measured from the ^^CO emission maps with no beam 
smearing. The opening angle appears slightly broader 
as the outflow axis becomes more parallel to the line-of- 
sight, but there is not a strong dependence on the in- 
clination. At low inclinations, the highest velocities are 
mainly perpendicular to the line-of-sight and the angle 
is measured using only a small sample of cells. When 
the outflow motion is mostly along the line-of-sight, the 
lobes become shorter and rounder. Observationally, out- 
flow orientations with completely overlapping lobes are 
not well suited to measuring opening angles compared 
to those with intermediate inclinations. Although the 
outflow inclination is difficult to infer from projection, 
the axes of the outflows investigated by AS06 are likely 
inclined between 30 to 60 degrees with respect to the 
line-of-sight. 

Figure [10] shows the opening angle versus time using 
the CO intensity maps to calculate the angle. Like Fig- 
ure m the angle increases on average as a function of 
time. Linear fits to the average opening angles in Figure 
Unigives = 0.80 + 0.18 log(t) and = 0.53 + 0.25 log(t) 
for inclinations of 30 and 45 degrees, respectively. How- 
ever, like Figure O the two lobes show different trends 
with time. For the projected data, the lower outflow lobe 
is nearly flat as a function of time, while the upper lobe 
widens significantly. This is mainly due to the concurrent 
widening and elongation of the lower lobe. Despite the 
variation, this suggests that trends inferred from molecu- 
lar ppv information are reflective of the actual ppp infor- 
mation. In this case, the similarity is particularly strong 
because CO is an efficient tracer of the low-density gas 
within the outflow cavity. 

4.4. Outflow Mass Evolution 

Molecular hydrogen has no dipole moment and is ob- 
servationally in visible unless the gas is strongly shocked 
(|Yu et al.lll999[ ). Consequently, CO, which is the next 




-2.0 -1.2 -0.5 0.3 1.1 1.8 2.6 

Fig. 8. — Log CO intensity (K) for Rl at t = 48kyr. Inclination 
with respect to the line-of-sight decreases vertically. The diameter 
of each map is 0.15 pc. The left column shows maps including 
all gas, while the right column shows integrated mass only for the 
outflow gas (l-u — i;*! > 2km s"-*^). On the right, the data has been 
convolved with a 4" beam assuming that the outflow is located at 
250 pc. 
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Fig. 9. — Opening angle versus inclination (see Figure 7) as- 
suming that the outflow is observed in -"^^CO at a distance of 250 
pc. The color scale indicates the age of the outflow in units of 10^ 
years. Symbols for the different lobes have been offset for clarity. 
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Fig. 10. — Rl opening angles as a function of time, where 
the angle is calculated using the integrated ■'^■^CO(l-O) intensity 
map including emission from channels with \v\ > 2 kms"-*^ and the 
beamsize is 1" . The thick and thin lines indicate 45 and 30 degree 
inclinations, respectively, relative to the line-of-sight. 

most abundant molecular species, serves as the primary 
tracer of molecular gas. In this section we consider two 
methods for estimating molecular gas mass from CO 
emission. 

For an assumed abundance of CO, a conversion factor 
can be used to relate CO intensity to total gas mass. 
This quantity, known as the X-factor, is defined as 



N{H2) 



(2) 



where W{^'^CO) is the ^^CO(l-O) integrated intensity 
and N{H2) is the column density of molecular hydro- 
gen. The abundance of CO depends upon the local 
gas density, temperature, and UV field, so that in re- 
ality, the X-factor varies widely over an indivi dual cloud 
(jPineda et al.ll2QQ8l : Iclover fc Mac Lq^|2Q11[ ) However, 



this complicated chemistry is usually reduced to an av- 
erage X-factor. As discussed in section 4.3, we adopt a 
typical mean CO abundance from observations. Here, we 
use the simulated outfiows to assess how accurately the 
standard X-factor and other more complicated methods 
recover outfiow mass. 

We adopt Xco = 1-8 x 102^cm-2K-^kms-\ which is 
the m ean value for molecular clouds in t he solar neighbor- 
hood (|Dame et al.ll2QQl[ ). iPineda et aP (2008) found that 
this factor overstimated the gas mass by a factor of 45% 
compared to the estimated extinction ma ss. However, 
the ty pical mean abundances inferred by iPineda et al .1 
([2008h are a factor of ~ 10-30 higher than the value we 
assume for our molecular line transfer calculation, where 
the difference is due to lower expected abundances for 
dense co re and outflow ^as versus lower density cloud 



g as fe.g..[Carolan et allffoOSl ). 



I Pineda et al.l ()2008l ) found that the X-factor is most 
reliable when the gas is both optically thin and Ay > 
4. While the emission from the outflow gas alone is 
marginally optically thin (with increasing optical depth 
as the inclination approaches 90 degrees and lower ve- 
locities), the integrated CO (1-0) emission through the 
core envelope is optically thick with r > 10. This means 
that the emission is saturated and Xco, which assumes 
a linear relationship between emission and gas mass, will 

underestimate the outflo w g:as ma ss. 

[Ballv et al.. (1999. ) and lYu et all (|1999[ ) present a more 
nuanced approach for obtaining the outflow gas mass 
fr om CO e mission. Their methods are similar to that 
of lArce GoodmanI (|200lf ). henceforth AGOl, which we 
follow here. All three techniques improve upon the es- 
timation above by correcting for optical depth effects 
and including more of the low- velocity outflow gas mass. 
We outline the AGOl proceedure below and post-process 
the simulations accordingly (see AGOl and references 
therein). 

While the CO (1-0) emission may be quite optically 
thick, it is possible to utilize the more optically thin 

CO (1-0) emission to calculate the true optical depth 
and correct the CO (1-0) mass estimate. In the flrst 
step of the procedure, we spatially average over all the 
^^00(1-0) and ^^00(1-0) spectra and derive the ratio of 
the two lines, i^i2/i3, as a function of velocity. The resul- 
tant parabolic line ratio shape indicates that the opacity 
is not constant with velocity as many studies assume. 

We then perform a polynomial flt of i^i2/i3('^), which 
is constrained to reach a minimum at the cloud velocity. 
Following AGOl, who limit the velocity range to exclude 
a secondary coincident cloud, we limit the fltted velocity 
range to within 1.5 kms~^ of the minimum. Using this 
flt, we extrapolate to the high- velocity wings, where ob- 
servationally ^^CO is too weak to be detected. (Even in 
the noiseless synthetic observations, the ^^CO emission 
is negligable at velocities above a few kms~^.) 

We next derive an effective "^^CO main beam tempera- 
ture," T^{x^y^v)^ to estimate the ^^CO opacity. Obser- 
vationally, if the CO (1-0) emission is greater than twice 
the rms noise, then T^{x^y^v) can be used directly. 
Otherwise, the ^^CO main beam temperature must be in- 
ferred from ^^CO: T^^ix,y,v) = T^^{x,y,v) / Ri2,iz{v^). 
We adopt a rms noise value of 0.06 K per 0.2 kms~^ 
channel, similar to that used by AGOl. 
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If CO (1-0) is indeed optically thin, then its opacity 
may be estimated: 



Ti3{x,y,v) 



-In 1 



To/[exp(To/Tex)-l]-0.87 



(3) 

where Tq = hu/k = 5.29 and Tex is the excitation tem- 
perature under the assumption that CO (1-0) is opti- 
cally thick: 



5.53 



ln[l + 5.53/(Tpeak + 0.82)]' 



(4) 



Here, Tpeak is the peak temperature for the -^^CO emis- 
sion. We find Tex ~10-12 K, in good agreement with the 
known simulation ambient gas temperature of 10 K. 

Finally, we derive the column density of ^^CO for each 
pixel: 



Nisix, y, v) = 2.42 x 10^^(Tex + 0.88) : 



Tis{x,y,v)dv 



l-exp(-To/Tex)' 
(5) 

where dv is the channel width in km s . The total mass 
per pixel of area, A, is given by 



M{x,y,v) = mu2Nu2{x,y,v)A, 



(7) 



where mu^ = 2.72mH is the mean molecular weight and 
the conversion between the ^^CO column density and 
the molecular hydrogen column density is assumed to be 
Nb_2 = 1.4 X lO^A^is. Note that the exact value depends 
upon the local X-factor. When modeling the HH300 
outfiow in the Taurus star forming region, AGOl adopt 
A^Hs = 7x lO^A^is, a g eneric value for dense gas in Taurus 
that was obtained bv iFrerking et all (p.982). 

An additional step is necessary to distinguish between 
cloud mass and outflow mass. Summing over the total 
area, we obtain M{v) as shown in Figure [TTl AGOl found 
that the peak, which includes the bulk of the core mass, 
was well fit by a Gaussian. We find that a Gaussian 
provides a satisfactory fit in the case of Rl, where there 
is only a single outfiow, while the R2 data, in which there 
is a second "hidden" outfiow, a Gaussian does not well 
describe all of the high density gas (see Figure [TT]). We 
then subtract this fit from M{v) and integrate over all 
velocities to get the outfiow mass, i.e., the shaded area 
shown on Figure [TTJ 

Figure [12] shows the simulated outfiow mass compared 
with that derived using the X-factor and the AGOl 
method. For R2, Figure [12] displays the the total mass 
ejected by both stars. Although the axes of the two R2 
outfiows are not aligned, projection, beam resolution, 
close proximity, and the low-velocities of the outfiows 
would prohibit individual identification of the two out- 
fiows observationally. 

In the figure, the "total ejected mass" is the total 
mass launched by the protostar according to the outfiow 

^ Alternatively, for optically thin gas the column density can be 
expressed as 

2.42 X 10i4(Tex + 0.88)T^3 i 



Ni3{x,y,v) 



[1 - exp(-To/Tex)] X [J(Tex) - J(Tbg)] ' 



(6) 



where J(T) = To/[exp(To/T) - 1] and Tbg = 2.73 K is the back- 
ground temperature. We find that this lowers the derived mass 
estimate by ~10%. 



AGOl Estimoted Moss 
Core Moss (Goussion Fit) 




AGOl Estimoted Moss 
Core Moss (Goussion Fit) 




Fig. 11. — Rl (top) and R2 (bottom) gas mass versus velocity 
derived using the AGOl method. A Gaussian fit over the velocity 
range —0.75 kms"-*^ < v < 0.75 kms"-*^ (dashed line) represents 
the non-outfiow material. The shaded region shows the blue-shifted 
and red-shifted outfiow mass. See AGOl Figure 7c for comparison. 



model. This is an essentially numerical quantity rep- 
resenting the total mass deposited in the zones around 
the protostar as described in Section 2. This quantity 
is distinct from what an observer would identify as the 
"outfiow mass," which is defined on the basis of velocity 
and thus also includes entrained gas. At any given time, 
the launched mass includes gas that goes on to comprise 
a low- velocity, less collimated component of the outfiow 
and gas that mixes with the ambient gas becoming in- 
distinguishable even with full information.^ Since the 
launched material is placed on the grid near the pro- 
tostar by design, this total does not include entrained 
envelope gas. Once the wind is deposited in the grid it is 
impossible to differentiate launched gas from entrained, 
accelerated gas. In principal, due to the addition of en- 
trained gas the total mass of high- velocity gas, i.e., the 

^ Some code methodologies, such as smooth particle hydrody- 
namics, permit the user to track individual fiuid elements, but 
ORION does not have this capability. 
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outflow mass, could exceed the mass numerically ejected 
by the protest ars. 

For both runs, the mass derived from the raw simula- 
tion data with |'^| > 2 kms~^ and outflow orientation of 
45 degrees is significantly less than the total ejected mass 
by a factor of 5-10. This reinforces the point that the 
majority of outfiow gas in an inclined outfiow is simply 
not observable even with a perfect method for converting 
between CO emission and total gas mass. However, the 
integrated ejected outfiow mass does not include gas that 
is entrained and accelerated by the outfiowing material, 
so in principle, the simulated and observed gas could be 
larger than the total ejected mass. 

The gas mass estimated from the CO emission us- 
ing the X-factor appears to track the mass of the high- 
velocity gas fairly well, although it is generally a factor 
of 5 — 10 lower. This discrepancy is due to the high mean 
optical depth of the core (r > 10 for CO (1-0) and r ~ 
a few for 1^00(1-0)). The AGOl method, which sub- 
tracts off the core mass, appears to perform somewhat 
better for most simulation outputs. However, the accu- 
racy of the mass estimation strongly depends upon how 
well M{v) can be fit by a simple function. For example, 
at late times the Rl M{v) is slightly asymmetric, leading 
to a negative derived outfiow mass. Note that Figure [TT] 
shows the log of the mass, whereas the fit is dominated 
by a few points around v km s~^. Differences be- 
tween the gas mass and fit at these low speeds far exceed 
the total mass contributed by high-velocity gas in the 
distribution wings. To account for this, we only subtract 
the curves beyond the full-width half-maximum of the 
Gaussian fit and omit velocities where the difference is 
negative. However, since the mass in the core is so much 
higher than the mass in the high- velocity channels, the 
mass estimates strongly refiect the symmetry and Gaus- 
sianity of the mass distribution. The sucess of AGOl in 
estimating the mass of outfiow HH300 relies upon the 
goodness of fit of a Gaussian and the relatively larger 
outfiow mass relative to the core mass. 

If the outfiow mass is instead defined as a sum over 
M{v) where \v\ > 2 kms~^ then the mass estimate fol- 
lows a similar trend to the mass derived using the X- 
factor. However, this suggests that the more compli- 
cated AGOl procedure still underestimates the total out- 
fiow mass in the optically thick case. The factor of ~5 
disrepancy between the two methods is due to different 
assumptions about CO abundance. 

4.5. Outflow Temperature 

Observationally, i2cO(l-0) and i2cO(2-l) may be 
used in combination to derive the excitation tem- 
perature, Tex- Temperature variation may be used 
to discrimin ate between different outfiow models 
(jArce fc Goo dman 2002). For example, shocked outfiow 
gas should be higher temperature than ambient gas, and 
the gas temp erature should r ise with the maximum out- 
fiow velocity (|Lee et al.ll200lf ). For optically thin gas, the 
excitation temperature may be related to the ratio of the 
line intensities by 



R 



'21/10 



4e 



-ll/Te, 



(8) 



where i^2 i/io is the ratio of the -^^C0(2-l) to 1^00(1- 
0) lines (jArce fc Goodmaj] l2002f ). Observers typically 
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Fig. 12. — Star and outflow gas mass as a function of time for 
the outflows inclined approximately 45 degrees with respect to the 
line-of-sight. The total ejected mass (red dotted line) includes all 
mass ejected over all velocities in the outflow. The outflow mass 
(green dot-dashed line) is estimated from the raw simulation data 
for \v\ > 2km s"-*^. The Xco mass (blue dashed line) is estimated 
from the ■'^^CO(l-O) emission observed with a 4" resolution beam 
in channels with \v\ > 2 kms"-*^. The thick purple dot dot-dot- 
dashed line shows the mass estimated via the AGOl method where 
a flt to the core gas mass is subtracted (see Figure [TT]) . The mass 
estimated from the AGOl method where a simple velocity cut of 2 
kms"-*^ is imposed in lieu of a flt is also shown (thin purple dot-dot- 
dot-dashed line). Note that the R2 mass estimates include both 
stars and both outflows. 
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assume that the high-velocity gas that constitutes the 
outflows is optically thin. We find that in many cases line 
ratios in the channels with \v\ > 2 kms~^ are above 1.0, 
the line ratio limit in the optically thick case, indicating 
that the gas is at most marginally optically thick in the 
outflow. 

Figure [13] shows that the estimated excitation temper- 
atures are generally higher than the ambient gas temper- 
ature of 10 K. The excitation temperatures also increase 
with time, which is expected if the shock strength and, 
hence post-shock temperature, increase with time. This 
effect also occurs in the simulations, although a rising gas 
temperature also results from increased protostellar heat- 
ing (see OKMK09). The two causes are difficult to distin- 
guish from a single average temperature value. However, 
Tex does not appear to be well correlated with the actual 
gas temperatures in the simulation. The disagreement at 
later times is likely due to the increas i ng op tical depth 
of the high- velocity gas. iHatchell et al.l (| 19991 ) show that 
gas with r = 1 and inferred Tex = 30 K under the as- 
sumption that the gas is optically thin will actually have 
Tex = 40K. The optically thin approximation increasely 
underestimates temperatures for higher optical depths 
and gas temperatures. 

5. DISCUSSION 

Often only a single outflow lo be is observed 
(jArce et al.l 120101 : iGinsburg et~aIll201lL creating uncer- 
tainty whether high gas velocites are instead due to co- 
herent turbulent motion. We find that even with com- 
pletely symmetric bipolar energy injection, interactions 
between the outflow and the turbulent envelope can re- 
sult in significant asymmetry between the lobes. Com- 
bined with inclination effects, it is therefore quite possi- 
ble for observations to identify only one component. 

Scale & Looney (2008), who include 12 outflows with 
only a single visible lobe in their photometrically iden- 
tified sample, attribute asymmetry to inclination effects 
that increase obscuration of the far lobe by the core en- 
velope. For outflows identified by high- velocity, coher- 
ent motion, the absence of a counterpart may be due 
to a combination of significant source advection, inter- 
ference with the turbulent cloud velocities, or asymmet- 
ric dynamical or magnetic interaction with the envelope. 
Since dense, star-forming gas and young protostars are 
observed to move subsonically relative to t heir en velopes 
and the host cloud (jKirk et al .1120101 : 1 Off ner et al.ll200 9a). 
obscuration and envelope interaction are the most likely 
causes of non-detected lobes. 

Outflow-envelope interactions are also largely respon- 
sible for the inferred outflow opening angle. Since the 
launching region and angle of the outflow remain fixed 
in our simulations, the significant angle broading is a re- 
sult of the outflow sweeping up sucessively larger solid 
angles of the core envelope. Over the course of the sim- 
ulation, the mean density of the core does not decrease 
significantly. However, the outflow mass and momentum, 
which are coupled to the protostellar mass, increase with 
time. Once the outflow breaks out of the dense envelope, 
an increasing amount of gas on the edge of the cavity is 
swept up, broadening the cavity. Thus, current obser- 
vational resolution is likely not yet probing the outflow 
launching region. Our hydrodynamic simulations are 
also in surprisingly good agreement with opening angles 



inferred from observations, suggesting that the details of 
the magnetic field evolution may play a sub-dominant 
role on these scales. Consequently, evolving opening an- 
gle trends can not yet be used to distinguish between 
different theoretical models and instead reflect a generic 
characteristic of outflow-envelope interaction. 

Exactly how outflows are responsible for driving turbu- 
lent motions in molecular clouds remains an open ques- 
tion in star formation. Some simulations of ~pc size 
clouds indicate that the momentum from outflows is suffi- 
cient to achieve quasi-steady state equip artition between 
turbulent and gravitat ional energy (|Nakamura fc Lil 
120071 : I Wang etM l2010f ). However, other simulations 
of jets indicate that wellcollimated outflows drive su- 
personic turbulence very inefficiently, although they 
may contribute sign ificant energy at lower velocities 
(jBanerjee et al.l 120071 ). However, observations of some 
molecular clouds suggest that the energy injected by out- 
flows is insufficent by an order of ma gnitude or more 
than the o bserved cloud turbulence (Arce et al." "2010"; 
IGinsburg et al. 20 U, but see also Swift & Welch 2008; 
Naka mura et al.ll201lf ). Such estimations include correc- 
tions for inclination effects and mass underestimations 
that amount to a factor of 4. Previous studies have found 
that neglecting low-velocity material may result in un- 
der predicting the outflow mass by at least a factor of 
2 (iMargulis fc Ladal[T985h . For outflows identified with 
minimum velocities of 10 kms~^, the missing momen- 
tum may be even larger since outflows from the lowest 
mass sources and outflows aligned perpendicular to the 
line-of-sight may be missed altogether. 

The method of AGOl, which includes lower velocity 
gas, exhibits the best agreement with the launched mass. 
However, it is unclear how much of the low velocity mate- 
rial, which dominates the mass estimate, is actually non- 
Gaussian turbulent core gas. The underestimation of 
outflow mass by both CO methods that employ a veloc- 
ity cutoff underscores the large uncertainties underlying 
these derivations. Such techniques are likely more accu- 
rate at later stages when the core mass has declined and 
more gas has been swept up by the outflows. Nonethe- 
less, this result has important implications for the deriva- 
tion of outflow momentum. In a worst case scenario, the 
masses of young, inclined outflows in optically thick re- 
gions are underestimated by factor of 50. 

Despite the severity of the problem, an order of magni- 
tude discrepancy relative to corrections already applied 
by observers is unlikely. In addition, the momentum 
contributed by outflows from very young low-mass stars 
is likely small compared to outflows from higher mass 
sources. Consequently, the discrepancy we demonstrate 
here is unlikely to account for the turbulent deficit mea- 
sured in clouds > 10 pc. In addition, outflows must 
not only account for the magnitude of th e energy injec- 
tion but th e appropriate range of scales. (Matznerl 120021 : 
ISwift fc W elch 2008). If measured outflow extents range 
from ~0.1-2 pc then they will not contribute driving mo- 
tions on larger scales. Magnetohydrodynamic simula- 
tions of outflow driven turbulence confirm that the peak 
of the velocity power spectrum occu rs at the maximurn . 
length scale of the driving outflows ([Carroll et al.ll2009l 
120101) . 

6. CONCLUSIONS 
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Fig. 13. — Outflow gas temperature versus time when the out- 
flow is inclined 45 degrees with respect to the hne-of-sight. The 
Rl and R2 gas temperatures are computed as an average over the 
cells with velocities \v45\ > 2km s"-*^. The CO excitation temper- 
ature, Tex, is calculated using equation [8] for an observation with 
4" resolution. 

We present results from self-gravitating, radiation- 
hydro dynamic simulations of low-mass protostellar out- 
flows using the adaptive mesh reflnement code, ORION. 
We produce synthetic observations of the simulations in 
^^CO and ^^CO in order to make direct comparisons with 
observations. 

We formulate a quantitative prescription for measuring 
outflow opening angles and test that it reproduces pre- 
vious estimates of observed outflow angles determined 
"by eye." Using this method, we flnd that both the 
simulation data and synthetic observations of the sim- 
ulation data show opening angle evolution that is similar 
to trends found by Arce & Sargent (2006). However, we 
show beam resolution can signiflcantly broaden the ob- 
served outflow angle a nd may broaden angles inferred by 
lArce fc Sa rgent^ ('2006') by as much as a factor of two. Dif- 
ferent inclinations with respect to the line-of-sight have 
a smaller effect on the measured angles. 

We flnd that the interaction between the outflows and 
the turbulent core envelopes produces signiflcant asym- 
metry and variation in outflow properties. Outflowing 
gas running into denser gas may be diverted or con- 
flned. This explains why observations may sometimes 
only identify a single outflow lobe. The clumpy outflow 
geometetry in some observed outflows may even result 
from outflow-envelope intereaction rather than variable 
mass accretion/ejection. 



Using CO isotopologues we estimate the observed gas 
masses in two ways. First, we infer mass using Xco, 
the standard factor often used by extragalactic observers 
to convert between CO emission and mol ecular gas 
mass. Second, we apply the method of Arce fc Good manI 
(|2QQT[ ). who combine ^^00(1-0) and ^^00(1-0) data to 
obtain mass estimates. We flnd that although the mass 
of the high- velocity outflow component tracks the actual 
outflow mass, this gas is not a good estimate for the 
total integrated outflow mass. Even with a perfect con- 
version between CO emission and gas mass, outflow mass 
inferred only using emission from the high velocity com- 
ponent underestimates the mass launched by a factor of 
- 2-5. 

Masses derived from the CO emission, which is a poor 
tracer of gas mass in the optically thick limit, underesti- 
mate the actual outflow mass by an add itional factor of 
5-10. Ideally, as I Arce fc GoodmanI (|200lf ) suggest, lower- 
velocity outflow gas can be included in the outflow esti- 
mate if the core mass can be modeled and subtracted 
from the emission. However, in the simulations the core 
mass distribution does not necesarily have a symmetric 
distribution of gas velocities, and so it is not well mod- 
eled with a simple functional form. Consequent l y, the 
more complicated method of 'Arce fc GoodmanI (|2001f ) 
produces a larger but not necessarily more accurate out- 
flow mass. 

We flnd that the excitation temperatures derived from 
^^CO are not necessarily well correlated with the simula- 
tion gas temperatures. However, they indicate that the 
high- velocity outflow gas is hotter than the ambient gas 
and that the gas temperature increases with time, which 
is consistent with the raw data from the simulations. 

In our simulations, the outflow launching region re- 
mains flxed throughout. Consquently, we conclude that 
observations of opening angle evolution do not probe 
the outflow launching mechanism and instead reflect the 
evolving interaction between the outflow and core enve- 
lope. Higher resolution observations in addition to sim- 
ulations including the effects of magnetic flelds will be 
needed to accurately investigate the physics of this in- 
nermost region. 
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